MODULE sefmm_cholesky_mo

   USE fmm_global_paras
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: sefmm_get_cholesky_MOs

   REAL(REALK), PARAMETER :: sefmm_cholesky_threshold = 1e-10

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE sefmm_get_cholesky_MOs(nBas,dens,CMO,rank)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: nBas
      REAL(REALK),   INTENT(IN)  :: dens(:)
      REAL(REALK),   INTENT(OUT) :: CMO(nBas,nBas)
      INTEGER(INTK), INTENT(OUT) :: rank

      REAL(REALK)   :: sq_dens(nBas,nBas)
      REAL(REALK)   :: Lij(nBas,nBas)
      REAL(REALK)   :: Pij(nBas,nBas)
      REAL(REALK)   :: tmp(2*nBas)
      REAL(REALK)   :: T0, T1
      REAL(REALK)   :: tol = sefmm_cholesky_threshold
      INTEGER(INTK) :: pivot(nBas)
      INTEGER(INTK) :: i,j, ij
      INTEGER(INTK) :: stat

      CALL CPU_TIME(T0)

      sq_dens = zero
      ij = 0
      DO j = 1, nBas
         DO i = 1, j
            ij = ij+1
            IF (ABS(dens(ij)) < sefmm_cholesky_threshold) CYCLE
            sq_dens(i,j) = dens(ij)
            sq_dens(j,i) = sq_dens(i,j)
         END DO
      END DO

      CALL fmm_matrix_norm('sq_dens from AO dens: ',sq_dens,nBas*nBas)

!      CALL DPOTRF('L',nBas,sq_dens,nBas,stat)
!      print *, 'DPOTRF status =', stat
!      CALL DGETRF(nBas,nBas,sq_dens,nBas,pivot,stat)
      CALL LEV3PCHOL('L',nBas,sq_dens,nBas,pivot,rank,tol,tmp,stat)
      print *, 'LEV3PCHOL rank =', rank
      IF (stat /= 0) THEN
         print *, stat
         !IF (stat > 0) print *, 'Dmat not positive definite!'
         CALL fmm_quit('sefmm_get_cholesky_MOs')
      END IF

      Lij = zero
      DO j = 1, rank
         DO i = j, nBas
            IF (ABS(sq_dens(i,j)) < sefmm_cholesky_threshold) CYCLE
            Lij(i,j) = sq_dens(i,j)
         END DO
      END DO

      Pij = zero
      DO i = 1, nbas
         Pij(i,pivot(i)) = one
      END DO

      CMO = zero
      CMO(:,1:rank) = IROOT2*(MATMUL(TRANSPOSE(Pij),Lij(:,1:rank)))

      CALL CPU_TIME(T1)
      CALL TIMTXT('>>> TIME USED get Cholesky MOs', T1-T0, LUPRI)

   END SUBROUTINE sefmm_get_cholesky_MOs

!-------------------------------------------------------------------------------

END MODULE sefmm_cholesky_mo

